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Abstract 

We present and discuss the development of an unconditionally stable algorithm 
used to solve the evolution equations of the Phase Field Crystal (PFC) model. This 
algorithm allows for an arbitrarily large algorithmic time step. As the basis for 
our analysis of the accuracy of this algorithm, we determine an effective time step 
in Fourier space. We then compare our calculations with a set of representative 
numerical results, and demonstrate that this algorithm is an effective approach for 
the study of the PFC models, yielding a time step effectively 180 times larger than 
the Euler algorithm for a representative set of material parameters. As the PFC 
model is just a simple example of a wide class of density functional theories, we 
expect this method will have wide applicability to modeling systems of considerable 
interest to the materials modeling communities. 
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1 Introduction 



The dynamics of a non-equilibrium system often results in highly complicated 
domain structures (microstructures). Typically, as time proceeds, the average 
size of these structures grows as a direct consequence of free-energy reduction: 
the interface is eliminated resulting in an increase in the size of homogeneous 
regions. Traditional non-equilibrium dynamics usually deals with the equi- 
librium states that are spatially uniform [Tf2ll3H] . i.e., the stable phases are 
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characterized by homogeneous values for the appropriate intensive thermody- 
namic variables. Classic, albeit quite simple, examples of models governing 
the evolution of such systems are the Cahn-Hilliard (CH) Equation for con- 
served systems [5] and Allen-Cahn (AC) Equation for non-conserved systems 
[6]. Examples are found in polymer mixtures [7], alloys [8.9J, liquid-crystals 
[TOfTT] . and in cosmology |12j . 

A model that has generated considerable recent interest is the Phase Field 
Crystal (PFC) Equation [TSUI], which is a conservative form of the famil- 
iar, non-conserved, Swift-Hohenberg (SH) Equation [13] . These systems differ 
from the CH and AC systems in that the stable phase is periodic. For SH 
models, the order parameter is viewed as capturing the inhomogeneities in 
a fluid associated with Rayleigh-Benard convection. In the case of the PFC 
model, which is a simple version of more elaborate density functional theories 
of liquid/crystal interfaces [ToTfT?] . the model captures features at the atomic 
scale, and thus contains highly detailed physical information about the sys- 
tem's structure. Such models can describe many of the basic properties of 
polycrystalline materials that are realized during non-equilibrium processing. 

The equations of motion governing these non-equilibrium phenomena are non- 
linear partial differential equations that cannot generally be solved analytically 
for random initial conditions. Therefore, computer simulations play an essen- 
tial role in our understanding and characterization of non-equilibrium phe- 
nomena. The standard Euler integration is known to be unstable for time step 
At above a threshold fixed by lattice spacing Ax [18]. In CH and AC systems, 
to maintain an interfacial profile, the lattice spacing must be smaller than the 
interfacial width £, and in PFC and SH systems, Ax must smaller than the 
periodicity selected by the system. Thus, the Euler update is inefficient, and 
in practice it is computationally costly to use to evolve large systems. Var- 
ious computational algorithms |19f20|21] have been developed by increasing 
At compared to the simplest Euler discretization. However, these methods 
still require a fixed time step, so they eventually become inefficient. Recently, 
unconditionally stable algorithms [22,23,24] were developed to overcome this 
difficulty for CH and AC Equations. These algorithms are a class of stable 
algorithms free of the fixed time step constraint for equations with a mix of 
implicit and explicit terms. While these algorithms allow for an increasing 
time step in CH systems as time proceeds, only a finite effective time step is 
possible for AC systems. A recent study [25J, based on this unconditionally 
stable algorithm, demonstrated analytically that one can use an accelerated 
algorithm At = At 2 ' 3 to drive the CH Equation, with the accuracy in corre- 
lation controlled by y/A. 

In this manuscript we apply this unconditionally stable algorithm to the PFC 
and SH Equations (Section 2). In Section 3 we establish the effectiveness of 
this approach through numerical studies of the algorithm, demonstrating that 



2 



the algorithm is both efficient and accurate for solving PFC Equation. Finally, 
in Section 4 we provide some concluding remarks. 



2 Unconditionally stable algorithms for PFC Equation 

In this section, we develop a class of unconditionally stable time stepping 
algorithms (At taken arbitrarily large without the solution becoming unstable) 
to the PFC and SH Equations. Although the main purpose of this section is 
to study unconditionally stable algorithms for the PFC Equation, we include 
a parallel discussion of the SH Equation, as the methodology applies to both 
equations with only trivial differences. 

2.1 Unconditionally Stable Finite Differences 

Both the PFC and SH Equations start from a free energy functional that 
describes the configurational cost of periodic phases in contact with isotropic 
phases, and can be expressed as 



where the periodic order parameter 0(x, t) has the wave number fc = 1 in 
equilibrium, and r < characterizes the quench depth. For the PFC equa- 
tion, r is proportional to the deviation of the temperature from the melting 
temperature T M — T. 

In the PFC model the order parameter (the density) is conserved, and thus 
the equation of motion is in the form of a continuity equation, d(p/dt = — V j, 
with current j = —MV(5F/5(p), where M is the mobility. Absorbing M into 
the time scale, we obtain the dimensionless form of the PFC Equation 



For the SH Equation, on the other hand, the order parameter is not conserved 
by the dynamics, and its evolution is postulated to have the the form 



F[*] = Jdx{\ 




(1) 



8_l = 5_F 

dt S4> 



V 2 {[r + (f + V 2 ) 2 ]0 + 3 }. 
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Eq. ([3]) has a simple dissipative form, where the rate of change of is propor- 
tional to the gradient (with an an L 2 inner product in functional space) of the 
free energy. 

In order to obtain an unconditionally stable algorithm, we now follow methods 
previously developed for the CH and AC Equations [2"31l2~4] . and work out in 
some detail how to semi-implicitly parameterize the equation of motion. We 
begin by "splitting" the linear terms in the equation of motion into "forward" 
and "backward" pieces, both for Eq. ([2]): 



<j) t+At + AtV 2 [(ai - l)(r + l)0 t+A * + 2(a 2 - l)V 2 t+At + (a 3 
<p t + AtV 2 k (r + l)0 t + 2a 2 V 2 0i + asV 4 t + 4 



H+At 

(4) 



and for Eq. (j3J) 



H+At 



-At 



[at - l)(r + l)<j)t+At + 2(a 2 - l)V 2 t+At + (a 3 - l)V 4 0* +Ai 



At 



ai(r + \)<k + 2a 2 V 2 t + a 3 V 4 ^ + 4>l 



(5) 



The constants a 1; a 2 and a 3 control the degree of splitting. In order to find 
the constraints on these parameters that yield an unconditionally stable algo- 
rithms, a standard von Neumann linear stability analysis on Eq. (J3]) and Eq. 
([5} may be performed. The procedures are quite similar and the results are 
identical for these two equations. We will only show the details for the PFC 
model in next subsection. 



2.2 Physical versus numerical instabilities 



As was found in the analysis of Vollmayr-Lee and Rutenberg [23] for the CH 
equation, the PFC equation will be linearly unstable to perturbations for le- 
gitimate physical reasons. Specifically, the isotropic phase 4> can be metastable 
or unstable to the stable periodic (crystalline) phase [14J if the system is an 
undercooled liquid. This situation (which is precisely what we are interested 
in modeling) is established when r + 3<ft 2 < 0. This physical instability com- 
plicates our standard von Neumann stability analysis, as we wish to predict 
when our numerical methods will cause an instability that is unrelated to the 
physical instability resulting from the thermodynamic. 

We can investigate the physical instability by a linear stability analysis on the 
equation of motion Eq. (|2"1). We let = <f) + i], where is a constant phase and 
7] is a small perturbation, and linearize the PFC Equation Eq. (|2l) in r\ to get 
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^ = v 2 



r + 



+ (1 + V 



2\2 



Vt- 



(6) 



This can be Fourier transformed to find 



Vk,t 



-k 2 



r + 



+ (l-k 



2\2 



(7) 



The physical instability for the above equation occurs for 



r k = k 2 (r + 30 2 ) + (1 - k 



2\2 



< o, 



which reduces to r + 3<p 2 < with k — 1 in the stable phase, as we indicated 
above. 

Now we can proceed to analyze the numerical stability and determine the 
constraints for the splitting parameters. We linearize the general step Eq. (j4j) 
by substituting (ft = + r\ and get 



7]t+At + AtV 2 (ai - l)(r + l)r/ 4+A i + 2(a 2 - 1)V 2 ^ +A4 + (a 3 - l)V% +A t 



= th + AtV 2 [ai(r + l)ry t + 2a 2 V 2 r/ t + a 3 V 4 r/t + 30 2 r/ 4 
The Fourier transform of the above equation results in 



(9) 



Vk,t+At 



1 - Atk 2 {( ai - l)(r + 1) - 2(a 2 - 1)A; 2 + (a 3 - l)k 4 } 



1 - Atk 2 { ai (r + 1) - 2a 2 A; 2 + a 3 fc 4 + 3cf> 2 } 



(10) 



This can be re-expressed as 



Note that r k = C k — 7?-k- While we want to avoid numerical instability, the 
physical instability is to be expected during the dynamics, and will not lead 
to numerical problems. But, as we indicated above, both of the instabilities 
will be captured by a general von Neumann stability analysis. One manner of 
dealing with this is to recognize that a proper unconditionally stable algorithm 
will be stable if and only if > and should be unstable if and only if < 0. 
The von Neumann stability criterion is |//k,t+At| < Vh.,t\- We can express our 
restriction on the regime of von Neumann stability as 



[1 + AtC k ] 2 > [1 + Atn k } 2 for r k > 
[1 + At£ k ] 2 < [1 + AtK k ] 2 for r k < 0. 
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The above inequalities can be rewritten as 



r k [2 + At(£ k + 1Z k )] > for r k > 
r k [2 + At(£ k + TZ k )) < for r k < 0, 

which, dividing by r k can be reduced to a single inequality, 2 + At(£ k + 7?. k ) > 
0, which implies < £ k + TZ^ for arbitrarily large At, and we obtain 



0<-k 2 



[r + l)(2ai - 1) + 30 2 - 2(2a 2 - l)k 2 + (2a 3 - l)k 4 , (13) 



which can be satisfied using the mode independent restrictions (and r > — 1) 



1 3<p 2 1 1 

ai < ; -, a 2 > -, a-i < -. (14) 

2 2(r + l)' 2' 3 ~ 2 v ; 

These are the constraints on the parameters a±, a 2 and 03 for unconditionally 
stable algorithms for all modes, for quenches in the range — 1 < r < — 30 2 . 
With these choices there is no threshold for At in order to maintain numerical 
stability. The quantity At is termed the algorithmic time step. We note that 
unconditional stability does not mean that the user of such algorithms may 
simply take as large a time step as is desired. Indeed, to obtain accurate 
physical results, there are additional restrictions on how large At may be. 



2.3 Effective time step 



To determine how large a time step we may take, and still maintain an ac- 
curate solution, we calculate the Fourier space "effective time step", as will 
be described below. We first note that when G l \ — (X2 — Q>3 — 1, Eq. (jU) 
corresponds to the traditional Euler update 



At Eu 



V 2 {[r + (1 + V 2 ) 2 ]^ + ^ 3 }, (15) 



where 4>' t+ ^ t denotes the field obtained after an Euler update on a previous 
field 4> t , while we use the unprimed <fit+At to denote the field obtained by 
unconditionally stable algorithm on cf) t throughout. 

We now define the spatial Fourier transform of kji = / dxe" ,k ' x 0<(x). In 
Fourier space, writing k 2 = |k| 2 , the Euler update becomes 



' } 'k,t+At - 0k,i 



At 



Eu 



fc 2 {[r+(l-* a ) a ]^ k|t +(0 8 ) k , t } l (16) 
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where (0 3 )k,t = / dxe * k ' x 0^(x). 



In Fourier space, the unconditionally stable algorithms Eq. can be written 
in a form that is analogous to Eq. ffTB"j) : 



<?k,t+Af — gkj 



-*»{ 



r + (1 - fc 



2\2 



c,t} , 



(17) 



where we define /c-dependent effective time step by 



AfPFC fu Af) = ^ 

e// 1 ' J l + AtF[(r + l)(l-ai) + 2A; 2 (a 2 -l) + A; 4 (l-a3) 
For SH Equation, the effective time step is 



At s J f (k, At) - i + At[(r + _ + 2p(fl2 _ 1} + fc4(i _ Q3)] . (19) 

At e ff(k, At) is an effective time step for a mode k, corresponding to an al- 
gorithmic time step At. Of particular interest in the case of periodic systems 
is the dominant mode (the lattice spacing in the PFC model), which, for the 
scaling choices made in Eq. (j2J) and Eq. ([3]) is simply k = 1. Using the param- 
eters employed in the simulations shown in the next section of r = —0.025, 
di = 0.45, ci2 = 0.5, a 3 = 0.5, we obtain the dominant effective time step for 
both equations 



At eff (k ,At) = — 

m °' ; l + 29At/800 

As At = oo, we obtain the maximum dominant effective time step At e ff(ko, oo) = 
800/29 ~ 27.6. We see that a large algorithmic time step At does not al- 
ways translate into a significant amount of system evolution, as the effective 
time step always remains less than 28 for these parameter choices, no mat- 
ter how large the algorithmic time step becomes. Thus, this value provides 
us with a useful bound on our exploration of just how large an algorith- 
mic time step to take, and still obtain accurate results. For example, if we 
take algorithmic time steps that yields an effective time step At e ff(k , At) = 
At eff (k , oo)/2 = 400/29, then we find At = 800/29. We demonstrate in the 
next section that this algorithm, when applied to the PFC Equation, real- 
izes a significant speedup compared to the traditional Euler algorithm, while 
maintaining a controlled level of accuracy. 



(20) 
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Fig. 1. Snapshots of simulations of the PFC model. Time increases from left to right. 
The first row shows the field obtained using the Euler algorithm with AtEu = 0.015. 
The second to bottom rows show the fields obtained employing the unconditionally 
stable algorithms, when using algorithmic time steps of At = 3, At = 10, and 
At = 30. 

3 Numerical results 

The simulations were performed in two-dimensions. Fig. [T] shows typical snap- 
shots of simulations for the PFC model with parameters = 0.07, Ax = 1.0, 
and L sys = 128 with random initial conditions which corresponds to the liquid 
state. For comparison, all the simulations start with the same initial condi- 
tion. In the Figure white regions indicate = 0, red = + 0.2 and blue 
= — 0.2. The top row was obtained using the Euler algorithm Atg^ = 0.015 
at time steps n = 30000, n = 60000, n = 90000, and n = 160000. The second 
and lower rows were obtained using the unconditionally stable algorithm with 
(moving down) At = 3, At = 10, and At = 30. For illustration and com- 
parison purposes, we show the system snapshots at the same energy density 
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Fig. 2. A measure of the error \f ((4> eu — 'Pun) 2 ) / ({4>eu — <fi) 2 ) versus the algorithmic 
time step At. 

as the top row — from left, the energy density E = 0.002374, E = 0.002360, 
E = 0.002357, and E = 0.002350 from the first to fourth column, respectively. 
We immediately see that, for the times and energies selected, there are no vis- 
ible differences between the Euler update simulation and the unconditionally 
stable algorithm with At = 3. However, there are visible differences between 
the Euler update and the simulations with At > 3. We now wish to make 
these qualitative observations more quantitative. 

To study the accuracy, we compare simulations at the same energy density 
E = 0.002374 (the first column in Fig. [1]). We compute a measure of the error: 



Y ((<fieu ~ 'Pun) 2 )/ {{4>eu — 4>) 2 ) , where e «( x ) denotes the fields obtained using 

Euler algorithm and </> un (x) denotes the fields obtained using the uncondi- 
tionally stable algorithm. Fig. [2] shows a plot of the error versus a range of 
algorithmic time steps At. Fig. [2] indicates that, unsurprisingly, the accuracy 
increases as we decrease the algorithmic time step At. When At < 3, the error 
is below 5%. On the other hand, the error behavior in Fig. [2] for a large algo- 
rithmic time step tends to saturate, mirroring the saturation in the effective 
time step At e ff for the dominant mode k = 1. 

Fig.[3]shows a comparison between the dominant effective time step At e ff(ko, At) 
in Eq. fTSOj) and a numerical estimate of the same quantity. The numerical es- 
timate is obtained by calculating t^*/n un , where tg* is the total time needed 
to reach the final state (a crystalline state without dislocations) using Euler 
algorithm and n un is the number of computer steps needed to reach the same 
state using unconditionally stable algorithms. We find good agreement for 
At < 3, while for At > 3, the separation between the analytic and numerical 
expressions increases. While the agreement at small times steps in unsurpris- 
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Fig. 3. A comparison between the theoretical dominant effective time step (solid 
line) and the numerical estimate of the same quantity (circle). 

ing, the curve provides a useful metric for the optimum algorithmic time step 
of At ~ 3, for the chosen parameters. When At = 3, the ratio of the number 
of time steps needed to achieve a particular energy using the unconditionally 
stable versus Euler algorithm is approximately 180 (the ratio of the dominant 
mode effective time step to the Euler time step). This is a substantial speedup, 
and requires minimal analysis to implement the technique. 



4 Conclusions 

In this paper, we have presented an unconditionally stable algorithm applica- 
ble to finite difference solutions of the the PFC Equation. We have demon- 
strated that a fixed algorithmic time step driving scheme may provide signif- 
icant speedup, with a controlled level of accuracy, when compared with Euler 
algorithm. For the representative parameters chosen, a speedup of a factor of 
180 was obtained. The analytical results and the numerical results are con- 
sistent with an effective time step analysis. Although this algorithm allows 
arbitrarily large algorithmic time steps, caution is indicated, as taking too 
large an algorithmic time step will yield inaccuracies with little improvement 
in the overall speedup of the calculation. This saturation in the speedup results 
from the details of how the system's energy evolution (and its corresponding 
microstructural evolution) is governed by the effective time step, which sat- 
urates as the algorithmic time step increases. Thus, there is little advantage 
in too large an algorithmic time step. A method for obtaining a reasonable 
value for the algorithmic time step At is suggested, in which a few test cases 
are run with different values of At to see which one offers a good speedup 
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and maintains the desired accuracy. The analytic form of the effective time 
step provides a useful guide for deciding how large a time step to select when 
trading off the obtainable speedup versus the loss of accuracy. 

We expect the methodology developed in this paper could find extensive ap- 
plications in a wide class of non-equilibrium systems. For example, it can be 
straightforwardly applied to the Swift-Hohenberg Equation, given its similar- 
ity with the Phase Field Crystal model. Additionally, the method also will 
apply to systems where there is a dominant mode realized at late times, such 
as is found in diblock co-polymers. This method should allow researchers to 
dramatically improve the computational efficiency associated with modeling 
the dynamics of materials systems. 
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